Numerical formulation of three-dimensional scattering problems for optical structures 
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This paper describes a numerical formulation for calculating wave propagation with high precision 
in a three-dimensional system. Yee's discretization scheme is used to formulate a frequency domain 
method that is compatible with the finite-difference time-domain (FDTD) procedure. When the 
S-matrix satisfies a unitarity (power flow conservation) condition, the method enables arbitrary S- 
matrix elements to be obtained within a numerical error of less than 10" 8 (2 x 10" 13 ) for double 
precision format. 

PACS numbers: 02.60.Cb,41.20.Jb,42.15.Eq,42.25.Bs 



I. INTRODUCTION 



Numerical simulations are important for studying 
electromagnetic-wave propagation in optical physics [l| 
and for designing silicon photonics (D, Q . One of the most 
successful numerical methods is the finite-difference time- 
domain (FDTD) method 0. It is suitable for visualizing 
dynamical propagation of electromagnetic waves. Simu- 
lations should not only give us a general understanding of 
the propagation, but also the details of the optical scat- 
tering. 

Designing chips such as silicon optical interposers [|J 
requires highly precise simulations including the trans- 
mittance and reflectance of the fundamental and higher 
order modes at each wavelength (i.e., each frequency). 
Small reflections may cause substantial instability [6| be- 
tween devices on the optical chip, and the small losses 
that result may build up (e.g. see Table I in Ref. [f|). 
Thus, we need a way to confirm that the numerical re- 
sults are precise when the scattering properties are sim- 
ulated. 

Here, we should note that the precision of a numerical 
calculation, which is affected by both numerical method 
(e.g. numerical stability for the FDTD |7|) and numeri- 
cal implementation (e.g. floating-point arithmetic |8|), is 
essentially different from the accuracy of numerical mod- 
eling that includes both the choice of the fundamental 
equation (e.g. microscopic nonlocal approach [9j is one 
such choice) and the discretization of the numerical pro- 
cedure (e.g. numerical dispersion for the FDTD @]). In 
a scattering simulation, the error related to the accuracy 
is often explicit and predictable, but the error related to 
the precision is apt to be implicit and unforeseeable. The 
S-matrix approach [lOj is widely used to study scatter- 
ing problems and numerical S- matrices have already 
been applied to scattering simulations on photonic crys- 
tal slabs [l2rll4j and metal films [la ]. Unfortunately, no 



method as yet has been discussed related to highly precise 
simulations in the three-dimensional optical structures. 

This paper proposes a numerical method that can pro- 
duce precise S-matrices for designing the silicon photon- 
ics devices. The method exploits a numerical procedure 
for quantum transport fl6L hj. The numerical precision 
of the method is evaluated in terms of the S-matrix prop- 
erties. 



II. FORMULATION 

Consider the macroscopic Maxwell equations in the an- 
gular frequency domain (u> space), 



V x H — —iLoe^e (x) E . 

V x E — zu^oM ( x ) H , 



(1) 



'Electronic address: 
URL: http://www.petra-jp.org/ 



t-usuki@petra-jp.org 



where Sq (/irj) is vacuum permittivity (permeability). 
The symbol i denotes an imaginary number, and the no- 
tation "exp (—iuit)" describes a harmonic oscillation. The 
symbols j, k, I, m, and n in the following formulation 
denote integers. Instead of using dipole moments in the 
optical media, Eqs. (Q]) are used to express the relative 
permittivity e (x) and relative permeability n(x). Note 
that Ime (x) , Im fj, (x) > for absorbing media. 

The optical system in Fig. Q] consists of three parts: 
two ideal waveguides and a region with scattering and 
absorption. The coordinates in Fig. [1] are ones trans- 
formed using x = (x, y, z) — > (u, v, w), and they are 
used to apply a non- uniform mesh to Eqs. ([1]). Note 
that the Jacobian matrix of this transformation is diag- 
onal in order to simplify the discussion in the following 
subsections A-F. 



A. Discrete representation 

Let us discretize the transformed space: (u, v, w) — > 
(I, m, n), and let us use cells in Yee's lattice [18j . Figure 
[2] shows an arrangement of discretized functions (e, /i, 
E and H) that are allocated a cell address (I, m, n). 
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where gk, and hk are denned as 



IB v(y) 



Figure 1: Optical system for multi-mode scattering. The 
propagation region is divided into three parts. 




Figure 2: A (I, m, n) cell in Yee's lattice. 



For example, the u component of the magnetic field 
H u (Z, m, n) means H u at u = Z + 1/2, v = m, w = n 
in the figure. The electromagnetic field in (x, y, z) coor- 
dinates is related to that of the (u, v, w) coordinates in 
the following manner. 



/i (0 ffo (m) h (n) H u (I, m, n) , 

fo (0 9i (m) ho (n) H v (I, m, n) , (2) 

fo (0 5o (m) hi (n) H w (I, to, n) , 



HoHy 

7m>h z 



and 



oE x 

EoEy 

£qE z 



fo (Z) gi (to) hi (n) E u (I, to, n) , 

fi (I) g (to) hi (n) E v (Z, to, n) , (3) 

fi (0 9i (m) Zz («) (') ™, ") j 



2 wdz 
/fc W _ c du 

Qu (m) = r- 

yfc v ; c cZw 



/\7 2 (n) 



c 



u=l+k/2 



v—m-\-k/2 



w—n+k/2 



(4) 



Here, c = l/^/eoMo is the velocity of light in a vacuum. 
We make the domain of variables u and v finite: < u < 
L and < v < M, where the boundaries L and M are 
integers. Furthermore, e v and /j, v for r\ — u, v, w in Fig. 
[2] are defined as 



£„ (Z, to, n) 
£« (Z, to, n) 
(Z, to, ra) 



lu=i, u=m+l/2, iu=n+l/2 ' 
\u=l + l/2,v=m,w=n+l/2 ' 
l«=;+l/2, «=m+l/2, w=n ' 



and 



Mm m,n) — [X \ u -i + i/2, v = m , w = n > 

fj, v (I, to, n) = ML=i !lJ=m+1 /2 )W =» . 



-m, w—n- 



-1/2 



The relative permittivity and relative permeability 
/i^ for i] = u, v, w satisfy periodic conditions, i.e., 
e n (Z + L, to, n) = e n (I, to, n) and e v (I, to + M, n) = 
(I, to, n). Figure [3] shows the coordinate transforma- 
tion (x, y) — > (u, v) described by Eqs. (|A3|) in Appendix 

m 

The six components of the electromagnetic field satisfy 
the following conditions: 

G(l + L, to, n) = B U Q (I, to, n) , 
G (I, m + M, n) = (Z, to, n) , 

where = H n . E v for rj — u, v, w. The parameters 
-B u , i?!, are generally complex numbers, and they satisfy 
I -B u | = \B V \ = 1. In this paper, B u = B v = 1. Now let 
us introduce the forward difference operators, 

A U G (I, to, n) =G(l+ 1, to, n) — G (Z, to, n) , 

A^^ (Z, m, n) = G (Z, to + 1, n) — ^ (Z, to, n) , (5) 

A„C/ (Z, m, n) = G (I, to, Ti+l) — (Z, to, n) . 

At it = L — 1 and u = M — 1, A u and A„ are defined as 

A u 5 (L - 1, to, n) = B U G (0, m, n) - G (L - 1, to, n) , 
A„g (Z, M - 1, n) = B„S (Z, 0, ») - S (Z, M - 1, n) . 

The backward difference operator is — Aj, where " T " de- 
notes the transpose. Using Eqs. (j4]) and ([5]), we can 
define modified difference operators: 

A u = /iA„/ , A„ = giA v g , A w = hiA w h . (6) 
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Here, 



H v = (H v (0, 0, n) ■ ■ ■ Hr, (L- 1, M - 1, n)) 1 , 
E r) = (£, (0, 0, n) ■ ■ ■ E V (L — 1, M - 1, n)) T . 

From Eqs. (Q]), the w-components H w in Eqs.© and E n 
in Eqs. ([3]) can be expressed as 



E w (n) = — ^- (A u i2"„ (n) - A V H U (n) ' 

£«, (ft) V y 

^ (ft) = -A-r (*lEv (n) - AJX (n)) . 
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Figure 3: Grid of Yee's lattice in periodic xy space: x (L) — 
as(0) = y(M) - y (0) = 4^m when L = M = 25, 
min (dx/du) = min(dy/dv) — 40 nm, and K u = K v — 5 (see 
Eq. (1A3[) in Appendix [A")l . (a) Grid with periodic boundary 
conditions (dot-broken lines), (b) Grid around the origin. 



B. Wave propagation in ideal waveguides 

For the bottom ideal waveguide, Mhe, M eh of Eqs. 
© and hj in Eqs. (gj are 

M HE (ft) = M bffi5 , M eh (ft) = M bBff , 
ft-o (n) = ft-i (n) = ft-h, as n < 0. 

For the top ideal waveguide, they are 

M HE (ft) = M tffB , M eh (ft) = M tjBH , 
/io (ft) — ^i (ft) = /i*, as 7i > iV — 1. 

Since the permittivity in the ideal waveguides is real, we 
have 



lmM KH E = lmM KE H = 



(9) 



for k = b, t. Figure g] depicts the arrangement of the 
above equations. The eigenvalue equation for the optical 



From Eqs. ([TJ) and ©, the u and v components of 
electromagnetic field satisfy 

A W H UV (n) = IMhe (ft) E vu (n) , 
-A^X U (") = iM £H (n) H uv (n) , 
with the 2LM x 2LM matrices, 



the 



(7) 



M„e = 



M 



EH 



-Aje~ 1 A„ H 
Aje^A, 



-Aje-iA, 



(8) 



and the LM x LAI diagonal matrices, 



e„ = diag ( e v (0, 0, n) , 
A*,, = diag ( Mi? (0, 0, ft) : 



.. , £,(1-1, M — 1, n) ' 
. . , Hr,{L- 1, M - 1, ra) 



The 2LM x 1 column vectors H uv and E vu in Eq. 
are expressed as 



)• 




Hum 



H v 



-E v 
E n , 



w(z) 

w = N — 

AM 



ra 
o 
in 



N-2 



o 
ra 



o 



Top waveguide channel 

^c'*'(0)...c«(ZJW-1) \ 
c}\0) . . . c,l'(LM-1) y J 2 



fc,(AM)=ft, EJW-1)imV w - 1 ) =0 
h (N-i) =h, HJ.NA) im e„(/v-1)=o 

h i( N -2) EJ/V-2) ^.,(W-2) 
h (W-2) H uv (W-2) ^(W-2) 
h,(N-3) EJW-3) ^,„(W-3) 



"o(2) H„„(2) 

h i( 1 > EJ1) 

"o( 1 > HJb 

h,(0) = h t E vu (0) 

*.(0) = »>, H uv (0) 



X c»(0)...c t '*>(/.M-1) \ 
c t »(0)...c t »(LM-1)y / T 

Bottom waveguide channel 




Figure 4: Graphic depiction of Eqs. ©-(O. 
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modes of the ideal waveguides is 

M kHE M kEH u k (j) = A 2 K (J) u K (j) , (10) 

for < j < LM. Here, u K are eigenvectors. Appendix 151 
show that all eigenvectors satisfy 

ul (j) M KEH u K (/) - 6 jjf = , (11) 

where is Kronecker's delta. The propagation con- 
stants (3 b and j3 t are given by 



Pk (j) = 2 arcsin 



. A s (i) 



2/7? 



for k = b, t , 



(12) 



and < arg/3 K (j) < n. We use an integer J K to sepa- 
rate the propagating modes and evanescent modes of the 
above j3 K : Im/3 K (j) = as < j < J K , and Im/3 K (j) 7^ 
as J K < j < LM. We build a square matrix consisting 
of the eigenmodes u K of Eq. (fTUf , a diagonal matrix 6 K 
from Eq. (fT2"|) consisting of the phases of the modes, and 

column vectors consisting of the coefficients c« (j) 
of the j-th mode in the waveguides: 

U K = (u K (0) • • ■ u K {LM - 1)) , 
B R =diag ( exp (i0 K (0)) , . . . , exp (ip K (LM - 1)) ) , 

^ = (<£=> (O)---cW (LM-1)) T . 

(13) 

Note that U^M kEH U k = 1 from Eq. (JTTJ). The (0) 
and (0) in the bottom ideal waveguide are defined by 
using Eq. (fT31 : 

(0) - 6 + 



E vu (0) = -^M bEH U b 



(14) 



1-0, 



(+) 



1-0, 



The H uv {N - 1) and (JV - 1) in the top ideal 
waveguide are defined in the same manner. 



C. Power flow with absorption 

Let us define the discretized formulation of time aver- 
aged power flowfl9j: 

Pz (n) = ^ Re (e\ u (n) fc x (n) ho (n) (n)) , (15) 

where "T" denotes the Hermitian conjugate. From Eqs. 
(HI]), (Unj) and (HU), the power flows P bz = P z (0) and 
P tz =P Z (N - 1) for Eq. HH) are given by 

Jk-1 



p** = E^U - 4 _) (i) 

J=0 



(16) 



7 K (i) 



4w 2 



cot 



0* U) 



where 7^ is real and positive. Equations ©, (|15p . and 
(UHJ) lead to the following relation between power flow 
and absorption: 

Pbz — Ptz 

3 N ~ 2 

= [ E ™ («) ( ImM ^ («)) ^™ (») 

n— 1 

+ flt„(n)(ImM M (n))H ue (n)l (17) 



c 

2^ 



JV-2 



(n) 

71=1 7]—U,V,W 

+ Hi (n) (Im^ (n)) H„ (r 
D. Transfer matrices 



We make 4LM x 1 column vectors V&fc of the electro- 
magnetic modes and 4LM x 4LM matrices Tk, which 
we call transfer matrices: *&k+i = T^k- The *f?k are 
defined as 



(+) 



E vu (n - 1) 



Equations ([7} and (|T4l yield the transfer matrices: 
( U b U b \ 



iM bEH U b 


1 iM bEH U b 1 




i-e b _1 ^ i-eo 





1 \ 


fto(n-l) 




/io(n) 


ho(n)/fi(n— 1) y 





1 \ 


/il(n-l) 


iM EH (n) , 1 < n 




hi(n)/i (n) / 



For T2N-1, we have 



-ih? (1 - e- 1 ) uj 

r - - 1 1 1 ^(i-or 1 )^ 



(18) 



(19) 



By using Eqs. (ITS)) and (IT9l . we can obtain the 2LM x 
2LM matrices t and f from the linear equation: 



)- T 2N-l---T ( I 



(20) 



which is the same as Eq. (2.17) in our previous study |16| . 
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E. Stable transfer matrix method 



F. Scattering matrix 



To solve the above Eq. (f20|) , we can establish a stable 
iterative procedure [H, [Tt} by using the 4LM x 4LM col- 
umn operator Pj . This procedure does not entail solving 
multi-slice eigenvalue problems for the region with scat- 
tering and absorption, which has advantages in both com- 
putational time and numerical precision over the other 
procedure [2(J that can be applied to optical scatter- 
ing [T^ - [l5| . In the following discussion, suppose we have 
2LM x 2LM blocks of matrices M and jV notated by 



M = 



The iterative equations for the 4LM x ALM matrix Ck 
and the 2LM x 4LM matrix Dk can be used to find t 
and f: 



Ck+i — TkCkPk, 

D k+ i = D k P k for < k < 2N - 1. 
The initial conditions are that 
1 



C " 1 1 

Ck always satisfies 

Cfe 



and D = ( 1 ) 



(21) 



(22) 



Cfc.oo Ckfli 
1 



because of the column operator Pk in Eq. (f2"Tj) . Pk has 
the following matrix representation: 



fe,10^fc,00 


1 



, , (23) 

T'fc,io^'fc,oi+T , fc,ii 

and here, the actual numerical procedure uses Gaussian 
elimination without partial pivoting for Pkii- From Eqs. 
and ([221), we find that iterating Eq. (gTJ) gives us 

i = C2NA0 and f = D 2 n, 00 • 

We can compute the electromagnetic field in the scatter- 
ing region by making 2LM x ALM matrices £ (n, k) for 
E vu and H (n, k) for H uv for the n-th cell. The initial 
conditions are 

£ (n, 2n + 1) = U (n, 2n) = ( 1 ) 

for 1 < n < N — 2. The £ (n, k) and H (n, k) are iterated 
using the column operator Pk'. 

£ (n, k + 1) = £ (n, k) P k for 2n + 1 < k < 2N - 1, 
H (n, k + 1) = U (n, k) P k for 2n < k < 2N - 1. 

Appendices [C] and |D] give two approaches to simulating 
the reverse scattering process from the top ideal waveg- 
uide to the bottom ideal waveguide. We can use either of 

these approaches to obtain i , r , £ (n, k) and li (n, k) 
for the reverse process. 



Let us define a J& X J& matrix r and a J t x Jt, matrix 
f only for propagating wave modes. The elements of r 
and t are normalized to f and t by the power flow of Eq. 
(fTo) as follows: 



7b (i) 

7b if) 



3 j 



t 



j j 



It (j) 
7b (j") 



jj 



If we only obtain the matrices r and t, we can reduce the 
matrix sizes of Ck and Dk'. 2LM x blocks Ck,oo and 
CVio, 2LAf x 2LM blocks Ck,oi and C^n, a, J b x J b 
block Dk. 10; and a J(, x 2LM block -Dfe,ii. By using 

the matrices f and t , we can define a J t x J t matrix 
r and a J& x Jt matrix t . In so doing, we obtain a 
(Jb + Jt) x (J b + J t ) S- matrix [ll|: 



S = 



r t 
t r 



(24) 



Let us define 2LM x J& matrices £ (n), H (n) 

(£00 {n,2N)^ 



(£(n)) jf 
(n (n)) ]f 



7b (j") 
(lioo (n, 27V) 



7b (j 



and 2LM x J t matrices £ , % only for propagating wave 
modes. We also define a 4LM x ( J& + J t ) matrix £ n and 
a ALM x 4LM matrix a„ 



a n = 



£(n) £ (n) 
•H (n) ft' (n) 7 ' 

c 3 / ImMjTB (n) 



2w 2 







ImMM(n) 



(25) 



From Eqs. (TTT1) and ([23)1 . the S-matrix including the case 
of absorption media satisfies 



N-2 



SlS-l+J2&<*nSn = 0- 



(26) 



This equation shows that the S-matrix is unitary when 
ImM HE (n) = ImM EH (71) = for < n < N - 1. 



III. NUMERICAL RESULTS 

The method described in the previous section is suit- 
able for analyzing very small scattering coefficients. Here, 
we will discuss the optical properties of a sidewall grat- 
ing waveguide (SGW) that is part of a phase shifter in 
a silicon optical modulator [22J. The grating structure is 
used to inject free carriers into the waveguide core, but 
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for it to work properly, the reflections of the fundamen- 
tal mode and radiation loss from the structure have to 
be suppressed. 

Figure [5] shows the SGW and its permittivity [23| dis- 
tribution. The SGW is a silicon waveguide and has a 
SiO*2 cladding layer on which a vacuum region is set. 
For the numerical analysis, we set the waveguide core 
to 440 nm x 220 nm and the grating pitch (width) to 284 
(74.5) nm. Before conducting the simulation, we have 




Figure 5: (a) SGW. (b) Permittivity distribution on a;?/-plane 
at z = 0. (c) Distribution on z:r-plane at y = 0. 

to obtain optical modes for the ideal waveguides at both 
ends of the SGW. The ideal waveguides have three waveg- 
uide modes and many radiation modes, as shown in Fig. 
[5J The propagating mode numbers and J t are that 
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Figure 6: Dispersion diagram of Si waveguide. There are 
three waveguide modes and many radiation modes in the SiC>2 
cladding layer and vacuum layer. 

Jb = J t = 125 — 96 at 1.4 — 1.6 fim, because we set four 
/im periodic boundary conditions along the x and y axes. 



The grid parameters of x and y are the same as in Fig. 
[31 whereas the grid parameters of z are dz/dw — 28.4 nm 
and N = 200. Figure [7| shows cross sections of the waveg- 
uide modes and radiation modes at a wavelength of 1.55 
/xm. The three waveguide modes in Fig. Eta)-(c) are la- 



<a^y=cl|pk^y^n§Fc^y^ 
■ HI - III o 



M "^ii*r^ ^111 y=ioo 



x(nm) 



xiiim) 



x (,um) 



Figure 7: Optical modes for Si waveguide at 1.55 /im. (a) j = 
0: Ef ± . (b) j = 1: E v xl . (c) j = 2: E^. (d) j = 3: radiation 
mode in the Si02 cladding layer, (e) j — 4: radiation mode 
in the Si02 cladding layer, (f) j — 100: radiation mode in 
the vacuum layer. 



beled Ef lt E\ x and E^i [19j- The waveguide modes are 
clearly localized around the silicon waveguide core. From 
the x-axis symmetry of the SGW, the fundamental mode 
Efi scattering through the SGW is intra-mode scattering 
(i.e. reflection) when the scattering only occurs between 
the waveguide modes. However, the scattering process 
between the Efy and radiation modes is complex. Fig- 
ure[7][d)-(f) shows three of the 100 radiation modes. The 
optical power outside the core is dominant for each radia- 
tion mode, but remains small inside the core. Therefore, 
we have to consider not only the reflection within the Ef x 
but also scattering between Ef t and other modes includ- 
ing radiation modes. 

Figure [5] shows the transmittance and reflectance of 
the SGW when Ef 1 is launched from the bottom. There 
is a stop band at around 1.45 /im, after which the 
fundamental-mode reflectance decreases with a periodic 
modification of increasing wavelength. The reflectance 
of each radiation mode also has a strong wavelength de- 
pendence. In particular, it has a low value at the point 
of the third mode reflectance in Fig. [5J The total opti- 
cal loss, which is caused by Sfj reflection and radiation- 
mode scattering, is 4.39% (— 13.6dB) at a wavelength of 
1.55 /im. Note that the loss does not have so strong a 
wavelength dependence. The simulation results clearly 
show that the sidewall grating does not cause significant 
losses or reflections for the silicon optical interposer Q. 

In an actual phase shifter [22], the gratings on either 
side of the waveguide core are doped with donors or ac- 
ceptors, and aluminum electrodes are attached to the 
gratings (see Fig. [5]). Doped silicon (25l. [261] with metal 
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Figure 8: Optical scattering from the fundamental mode. The 
SGW has a stop band at 1.45 /tm. This is a typical filtering 
property for |roo| 2 , and it is designed for a modulator with low 
optical loss in the range of 1531 — 1591 nm [24| . The radiation 
and reflection loss is about — 15 dB respectively. 



-type waveguide 



p -doped grating 

(B, 10 19 cm" 3 ) 



n -doped grating 
(P, 10 19 cm" 3 ) 




Al electrode 
(anode) 



Al electrode 
(cathode) 



Figure 9: P-i-n diode in SGW. The left (right) side of the 
side wall grating is doped with a donor (acceptor) density 
10 19 cm -3 . The electrodes' gap is 2.44 /im. The permittivity 
of doped silicon is taken from p5j. [2r3|. and the permittivity 
of aluminum is taken from [271. |28|. 



electrodes [27I . l28fl changes the optical index and absorp- 
tion properties of the SGW. We calculated the scatter- 
ing process for absorption at a wavelength of 1.55^tm. 
Figure [10] shows the distribution of |tjo| and |tjo| as 
< j < 103. The total optical loss is 5.59 % (-12.5 dB), 
and it consists of Ef x reflection (see the case of j = 
in Fig. [l"0|) . E21 scattering {j — 2), radiation-mode 
scattering (3 < j < 103), and absorption loss (0.52% 
(-22.8 dB)). Note that the distribution around -80 dB 
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Figure 10: Distribution of scattering coefficients on the SGW 
at 1.55/im. The distribution at around — 30 dB of the SGW 
with the p-i-n diode is not so different from that of the SGW 
without the diode. 



and the E^ scattering are caused by an x-axis asymme- 
try due to the different optical indexes of the p-doped 
and n-doped regions (2ol l26j. 



IV. CONCLUSIONS 

This proposed method can calculate scattering coeffi- 
cients for all incident modes from the bottom waveguide 
and the top waveguide. To determine the precision nu- 
merically, we use the max norm S max on the left side of 
Eq. ([261): 



W-2 



(27) 



and the max norm U K max on the left side of Eq. (JTTJ for 
the propagation modes: 



U K 



[u K (0)---«„ {J K -l)) f M KEH 
x (« B (0)».«„(J„-1))-1|I 



(28) 



We can estimate the numerical error of the obtained scat- 
tering coefficients by performing a double precision cal- 
culation of Eqs. (J27J an d ([25]). Figure [TT] shows that 
numerical error of the eigenvalue calculation for ideal 
waveguide modes causes S max error. Using the Type I 
reverse propagation described in Appendix [C] results in 
an S'max value of less than 10~ 8 . On the other hand, the 
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Figure 11: U K ma x dependence of 5 max - Circles (crosses) in- 
dicate Type I (II) results. Ub max = Ut max in this case. The 
inset plots the absolute value of the 00-element on the left 
side of Eq. flTTJ. 



numerical error of the power-flow conservation for Ef x 
scattering does not depend on U K max , and it is less than 
2 x 10 -13 , as shown in the inset of Fig. [TTJ Therefore, 
our method satisfies the condition placed on the S-matrix 
(Eq. |27p) and gives us detailed optical properties with 
high enough precision for designing silicon photonics de- 
vices [5. 22|. For typical simulations, we recommend the 
Type II in Appendix Q2 because it takes only about half 
of the computational effort of Type I. 

We will use this method to analyze low and complex 
optical scattering cases. Furthermore, the discretiza- 
tion of the permittivity distribution on the Yee lattice is 
compatible with FDTD. By combining this method and 
FDTD, numerical analyses with an optical propagation 
model can be made general and detailed. 
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Appendix A: Non-uniform mesh for equations ([4]) 

Here, we show the details of the coordinate transfor- 
mation. The periodic boundary conditions should be ap- 
plied to dx/du and dy/dv: 

dx , . dx , T . dy , . dy , 

— (u ) = —{u + L , / (v = JL (v + M) . 

du du dv dv 

We introduce the periodic function F (£, K) into dx/du 
and dy/dv. 



K-l 



F(£, tf) = 2 2 *-W*«) ]J 



J 



J=i 



K+J 



(Al) 



The function (|A1[) has several properties: 



F(£ + 1,K) - F K) = F (£, K) , 

K-l 

maxF(£, K) = F (0, K) = 2 2I< - 1 J] 



J 



j=i 



K + J' 



min F (C, K) 

[Fit, K)dt 
Jo 



hm F(t;, K) 

K — >cx) 



Fl-,K 



= 1, 



Z=-c 



where S (^) is the Dirac delta function. From Eq. 
(|Alj) . we obtain an analytical formula for the integral 
ofF(£, K). 



F 



•/=! \J'=1 



K + J' 



2irJ 



The integral of F has a staircase shape. For example, we 
can set x (it) and y (v) by using 



x (it) — x (0) — u min (dx/du) ru/L 



x (L) — x (0) — L min (dx/du) J 
y (v) - y (0)-v mm (dy/dv) f v/M 
y(M)-y(0)-M min (dy/dv) ~ J Q 



F (L K u ) dt 
(A3) 



given ten parameters L, M, x(0), 2/(0), x(L), y (M) 
min (dx/du), mm(dy/dv), K u and K v . 
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Appendix B: Orthogonality of eigenvalue equation 

(unj 



for N > n> 1 , 



We should note that u\ (j) M k eh u k (j) — 0, where 
"T" denotes the Hermitian conjugate, when ImA^ (j) ^ 0, 
and that no linear combination aM k eh + bM^^ (real 
a, b) is positive definite. However, the eigenvector u K 
still satisfies 

(j) M kEH u k (j) 0. 
As an example, let us consider a simple equation, 

is not positive definite. This equation has complex eigen- 
values and eigenvectors: X± = ±i, u± = ( } . | , and 



±i 



1 

-1 



u±^0. 



Therefore, we can normalize all eigenvectors as 

ul (j) M kEH u k (j') = S jjt 

when eigenvalues are non-zero and non-degenerate, and 
number of eigenvalues is equal to the order of M k eh- 
The eigenvalue equations for the numerical results of Sec. 
IIIII satisfy this condition. 

We recommend checking which eigenvalue equations 
satisfy the condition or not, because there exists a 

counter example: ^ ^ ^ M = ^ ( ^1 ^ U ' ^ 1 * s 
equation has only one linearly independent eigenvector 
1 



-1 



and 



J- 2n ~ I hi(n) -iM EH (n) I , 

\ M«-i) m«-i)M«) / 
f 1 \ 

T 2n -1 — M») -iM HE (n-\) 

\ h (n-l) h (n-V)hi{n-l) I 



T _ , -ih 2 b (1 - e b ) uj 



ihiu b (i-e b )u£ i 



The reverse equation corresponding to the forward Eq. 

(na is 



- T o-" T ^i(f.' 



Reverse iteration is defined as 

D k = D k+1 P k for 0<fc<2iV-l, 
with initial conditions, 

C ™ = ( 1 ) and D 2N = ( 1 ) . 

The column operator P k can be defined in the same man- 
ner as Eq. (j2"5)) ; that is, 



Z> ^ ^fc.lO^fc+l.OO 



-^ , fe.io^fc+i ] oi"f"-^'&.ii 



and 





l 



u 



2 1 
1 



u = 0. 



Accordingly, we should carefully normalize eigenvectors 
when the equation does not satisfy the condition. 



Appendix C: Reverse propagation type I 

For reverse propagation in the same framework as that 
of forward propagation, we define a ALM x 1 column 
vector *& k as 



1 

1 



*t for 2N >k>0. 



The ALM x ALM transfer matrix T k , which satisfies 
— T k ^ k+1 , is defined as 



2N-1 



iM tEH Ut 1 iMtEHUt 1_ 

i-a; 



hi i- 



without k = 2N —1,1. P 2 n-i an< ^ -f i are 

1 



P2N-1 — P\ 



1 



because T 2N _ 1 (T ) is different from To (T2jv-i)- Here, 

( u t u t \ 

= -iMtEHUt 1 -iMtEHUt 1 

\ hi(N-2)h t hi(N-2)ht l-O t J 



Thus. 



and 



(U t {l + 9 t ) ih 1 (N-2)h t U t (l-0 t )U 1 t 
1 
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Furthermore, 



Thus, 

t' c\ 



C, = T, ( C ^ 00 C 2,oi | /> 



1 



o -ihf (i - e b ) uj 

C'1,10 C' liU + i/igE7 6 (l-0 6 )tf2 



From Eqs. (fTO)) . T C 1 is an expression similar to 
T 2 n-iC 2 n-i- Reverse iteration yields 



At the 27V — 2 — > 27V — 1 step, we introduce C 2N _ 1 and 

•^2AT-l : 



'2N-1 



-IVi + C 2 JV-1,01 — TT 



i-e t 





C^iV-l.Ol 
1 



2JV-1 



^ n iM tEH Ut i n A 

^£>2JV-l,01 ff-^J^ D2N-1,01 J 



Thus. 



* — Cq.oo and r — -^o.io • 



At the n-th cell, we can introduce £ (n, 27V — k) and 
H' (n, 27V -k). For 1 < n < TV -2, the initial conditions 
are 

£ (n, 27V - 2n - 2) = %' (n, 2N - 2n - 1) = ( 1 ) . 

The following iterations derive £ (n, 2N) and 
# (n, 27V): 

£ (n, 2N-k) = £ (n, 27V - ft - 1) 

for 2n + 1 > ft > 0, 

(n, 2N-k) = H (n, 2N - k - 1) P k 
for 2n > ft > 0. 



Appendix D: Reverse propagation type II 

We can derive another equation for f and t : 

( \ 



C 2 N — T2N-lC 2N _ 1 P 2 N-li 



2N-1 



2N-2 ■ ■ ■ J- 



iMtEHUt 1 
hf l-0t 



(Dl) 



,D 2 N — D 27V-l-f , 2W-l- 

Note that C k m = C k .oi, D k,oi = D k,0i and P k 01 = 
Pfc.oi for 27V — 1 < k < 27V. The iterative procedure of 
Eq.' dDTJ) yields 



r = C 2N00 and t = D 



2N, 00 • 



We can add £ (n, 27V- 1) and H (n, 27V- 1) to the 
27V - 2 -)• 27V - 1 step. 



£(n, 27V -1) = (0 27V - 1) ) , 

H (n, 27V - 1) = ( -Hoi (n, 27V - 1) ) . 

Finally, we obtain 

£ (n, 27V) (n, 27V - 1) P 2N _ V 
H (n, 27V) = % (n, 27V - 1) P 2N _ V 
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